library(tidyverse)
library(RColorBrewer)
library(readbitmap)
library(rjson)
library(patchwork)
geom_spatial <- function(mapping = NULL,
data = NULL,
stat = "identity",
position = "identity",
na.rm = FALSE,
show.legend = NA,
inherit.aes = FALSE,
...) {
GeomCustom <- ggproto(
"GeomCustom",
Geom,
setup_data = function(self, data, params) {
data <- ggproto_parent(Geom, self)$setup_data(data, params)
data
},
draw_group = function(data, panel_scales, coord) {
vp <- grid::viewport(x=data$x, y=data$y)
g <- grid::editGrob(data$grob[[1]], vp=vp)
ggplot2:::ggname("geom_spatial", g)
},
required_aes = c("grob","x","y")
)
layer(
geom = GeomCustom,
mapping = mapping,
data = data,
stat = stat,
position = position,
show.legend = show.legend,
inherit.aes = inherit.aes,
params = list(na.rm = na.rm, ...)
)
}
This example shows one how to estimate immunofluorescence intensity at each spot. It is for educational purposes only and is not supported by 10x genomics.
We will be using a public data set from: https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Adult_Mouse_Brain_Coronal_Section_1
Note: all paths are relative to my home directory and will not work on your system. Adjust them as needed.
The image channel setup. This is different from what is on the support site. The support documentation is correct for the loupe file
tmpdir <- tempdir()
url <- 'https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Adult_Mouse_Brain_Coronal_Section_1/V1_Adult_Mouse_Brain_Coronal_Section_1_spatial.tar.gz'
file <- basename(url)
download.file(url, file)
untar(file, exdir = tmpdir)
list.files(tmpdir, recursive = TRUE)
## [1] "spatial/aligned_fiducials.jpg" "spatial/detected_tissue_image.jpg"
## [3] "spatial/scalefactors_json.json" "spatial/tissue_hires_image.png"
## [5] "spatial/tissue_lowres_image.png" "spatial/tissue_positions_list.csv"
images_cl <- read.bitmap("spatial/tissue_lowres_image.png")
images_chnl_1 <- images_cl[,,1]
images_chnl_2 <- images_cl[,,2]
images_chnl_3 <- images_cl[,,3]
Take a quick look at the image channels
image(images_chnl_1)
image(images_chnl_2)
image(images_chnl_3)
Get the height and width of the images
height <- data.frame(height = nrow(images_cl))
width <- data.frame(width = ncol(images_cl))
Make the images into grobs for easy plotting in ggplot2
grobs <- list(grid::rasterGrob(images_cl, width=unit(1,"npc"), height=unit(1,"npc")))
grobs_chnl_1 <- list(grid::rasterGrob(images_chnl_1, width=unit(1,"npc"), height=unit(1,"npc")))
grobs_chnl_2 <- list(grid::rasterGrob(images_chnl_2, width=unit(1,"npc"), height=unit(1,"npc")))
grobs_chnl_3 <- list(grid::rasterGrob(images_chnl_3, width=unit(1,"npc"), height=unit(1,"npc")))
images_tibble <- tibble(grob=grobs,
grob_chnl_1 = grobs_chnl_1,
grob_chnl_2 = grobs_chnl_2,
grob_chnl_3 = grobs_chnl_3)
images_tibble$height <- height$height
images_tibble$width <- width$width
images_tibble
## # A tibble: 1 x 6
## grob grob_chnl_1 grob_chnl_2 grob_chnl_3 height width
## <list> <list> <list> <list> <int> <int>
## 1 <rastrgrb> <rastrgrb> <rastrgrb> <rastrgrb> 600 600
Get the scale factors so that we can adjust the spot pixel positions for the lowres image
scales <- fromJSON(file = "spatial/scalefactors_json.json")
Make a merged dataframe with all the information we need.
bcs <- read.csv("spatial/tissue_positions_list.csv",
col.names=c("barcode","tissue","row","col","imagerow","imagecol"), header = F)
bcs$imagerow_scaled <- bcs$imagerow * scales$tissue_lowres_scalef # scale tissue coordinates for lowres image
bcs$imagecol_scaled <- bcs$imagecol * scales$tissue_lowres_scalef
bcs$imagerow_scaled_round <- round(bcs$imagerow * scales$tissue_lowres_scalef) # Rounded scales
bcs$imagecol_scaled_round <- round(bcs$imagecol * scales$tissue_lowres_scalef)
bcs$tissue <- as.factor(bcs$tissue)
bcs$height <- height$height
bcs$width <- width$width
Function to get intensities
This is taking into account the spot diameter and pulling the average of pixels in every direction Might be better if using the full resolution image but seems to perform well on the lower resolution images.
get_intensity <- function(bcs, image, scales) {
radius <- round((scales$spot_diameter_fullres * scales$tissue_lowres_scalef)/2)
intensity <- list()
for (i in 1:nrow(bcs)) {
temp <- expand.grid(x = (bcs$imagerow_scaled_round[i] - radius):(bcs$imagerow_scaled_round[i] + radius),
y = (bcs$imagecol_scaled_round[i] - radius):(bcs$imagecol_scaled_round[i] + radius))
intensity[[i]] <- mean(image[temp$x, temp$y])
}
return(intensity)
}
Calculate intensities for each channel and each spot
intensity_1 <- get_intensity(bcs, images_chnl_1, scales)
intensity_2 <- get_intensity(bcs, images_chnl_2, scales)
intensity_3 <- get_intensity(bcs, images_chnl_3, scales)
Add the intensity info to bcs
bcs$intensity_1 <- unlist(intensity_1)
bcs$intensity_2 <- unlist(intensity_2)
bcs$intensity_3 <- unlist(intensity_3)
head(bcs)
## barcode tissue row col imagerow imagecol imagerow_scaled
## 1 ACGCCTGACACGCGCT-1 0 0 0 2238 3201 55.39604
## 2 TACCGATCCAACACTT-1 0 1 1 2478 3337 61.33663
## 3 ATTAAAGCGGACGAGC-1 0 0 2 2240 3476 55.44554
## 4 GATAAGGGACGATTAG-1 0 1 3 2480 3612 61.38614
## 5 GTGCAAATCACCAATA-1 0 0 4 2242 3751 55.49505
## 6 TGTTGGCTGGCGGAAG-1 0 1 5 2482 3887 61.43564
## imagecol_scaled imagerow_scaled_round imagecol_scaled_round height width
## 1 79.23267 55 79 600 600
## 2 82.59901 61 83 600 600
## 3 86.03960 55 86 600 600
## 4 89.40594 61 89 600 600
## 5 92.84653 55 93 600 600
## 6 96.21287 61 96 600 600
## intensity_1 intensity_2 intensity_3
## 1 0.02352941 0.01176471 0.04329412
## 2 0.02352941 0.01176471 0.04313725
## 3 0.02352941 0.01176471 0.04454902
## 4 0.02352941 0.01176471 0.04313725
## 5 0.02352941 0.01176471 0.04313725
## 6 0.02352941 0.01176471 0.04313725
Define our color palette for plotting.
myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
I like this color scale but it’s easily changed with code like this
scale_fill_gradient(low = "black", high = "red")+
instead of
scale_fill_gradientn(colours = myPalette(100))+
Plot the data with and without spot intensities
plot1 <- bcs %>%
filter(tissue =="1") %>%
ggplot(aes(x=imagecol_scaled,y=imagerow_scaled,fill=intensity_1)) +
geom_spatial(data=images_tibble[1,], aes(grob=grob_chnl_1), x=0.5, y=0.5)+
geom_point(shape = 21, colour = "black", size = 2, stroke = 0.1)+
coord_cartesian(expand=FALSE)+
scale_fill_gradientn(colours = myPalette(100))+
xlim(0,max(bcs %>%
dplyr::select(width)))+
ylim(max(bcs %>%
dplyr::select(height)),0)+
xlab("") +
ylab("") +
labs(fill = "NeuN Intensity")+
ggtitle("NeuN Spot Intensity") +
theme_set(theme_bw(base_size = 10))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_blank())
plot2 <- bcs %>%
filter(tissue =="1") %>%
ggplot(aes(x=imagecol_scaled,y=imagerow_scaled,fill=intensity_1)) +
geom_spatial(data=images_tibble[1,], aes(grob=grob_chnl_1), x=0.5, y=0.5)+
coord_cartesian(expand=FALSE)+
xlim(0,max(bcs %>%
dplyr::select(width)))+
ylim(max(bcs %>%
dplyr::select(height)),0)+
xlab("") +
ylab("") +
ggtitle("NeuN Image Only") +
theme_set(theme_bw(base_size = 10))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_blank())
plot1 + plot2
plot1 <- bcs %>%
filter(tissue =="1") %>%
ggplot(aes(x=imagecol_scaled,y=imagerow_scaled,fill=intensity_2)) +
geom_spatial(data=images_tibble[1,], aes(grob=grob_chnl_2), x=0.5, y=0.5)+
geom_point(shape = 21, colour = "black", size = 2, stroke = 0.1)+
coord_cartesian(expand=FALSE)+
scale_fill_gradientn(colours = myPalette(100))+
xlim(0,max(bcs %>%
dplyr::select(width)))+
ylim(max(bcs %>%
dplyr::select(height)),0)+
xlab("") +
ylab("") +
labs(fill = "DAPI Intensity")+
ggtitle("DAPI Spot Intensity") +
theme_set(theme_bw(base_size = 10))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_blank())
plot2 <- bcs %>%
filter(tissue =="1") %>%
ggplot(aes(x=imagecol_scaled,y=imagerow_scaled,fill=intensity_2)) +
geom_spatial(data=images_tibble[1,], aes(grob=grob_chnl_2), x=0.5, y=0.5)+
coord_cartesian(expand=FALSE)+
xlim(0,max(bcs %>%
dplyr::select(width)))+
ylim(max(bcs %>%
dplyr::select(height)),0)+
xlab("") +
ylab("") +
ggtitle("DAPI Image Only") +
theme_set(theme_bw(base_size = 10))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_blank())
plot1 + plot2